Posterior sampling
========================

.. index:: rsamplers.rsample.sample

The sampling function
-----------------------

.. include:: rsamplers_sample.rst

.. index:: rsamplers.rwmh

Algorithm : Random-walk Metropolis Hastings
---------------------------------------------

.. include:: rsamplers_rwmh.rst

.. index:: rsamplers.imh

Algorithm : Independent Metropolis Hastings
------------------------------------------------

.. include:: rsamplers_imh.rst

.. index:: rsamplers.apt

Algorithm : Adaptive parallel tempering
------------------------------------------

.. include:: rsamplers_apt.rst

.. index:: rsamplers.slice

Algorithm : Slice Sampler
------------------------------

.. include:: rsamplers_slice.rst

.. index:: rsamplers.usrsmplr

Algorithm : usrsmplr Sampler
------------------------------

.. include:: rsamplers_usrsmplr.rst

User-defined Algorithms: Another approach
-------------------------------------------

Due to the fact that the RISE toolbox is written in a modular manner, a user can apply his own sampling algorithm. However, if he wants RISE to
process the draws, the output of the sampling should be organized in a way
that RISE understands. Here is an example in which we write a wrapping function **wrapped\_dramrun** around
the delayed-rejection adaptive Metropolis (DRAM) sampler available `here <https://mjlaine.github.io/mcmcstat/>`_ .

.. code-block:: matlab
   :linenos:

	function r=wrapped_dramrun(m,x0,SIG,N)
	% wrapped_dramrun : RISE wrapper for the DRAM code
	
	[~,lb,ub,x0_,vcov]=pull_objective(m);
	
	[model,data,param,opts]=dramrun_inputs();
	
	% call to the DRAM code
	[results, chain] = dramrun(model, data, param, opts);
	
	r=results2rise();
	
	    function [model,data,param,opts]=dramrun_inputs()
	        
	        if isempty(SIG),SIG=vcov; end
	        
	        if isempty(SIG),x0=x0_; end
	        
	        % RISE assumes column vector, DRAM assumes row vector
	        param = struct(); param.par0 = x0.';  param.bounds=[lb,ub].';
	        
	        opts=struct(); opts.qcov = SIG; opts.nsimu=N;% Number of samples
	        
	        data = {};
	        
	        model=struct();
	        model.ssfun = @(p)m.routines.likelihood(p.',m);
	        model.priorfun = @(p)log_prior_density(m,p.');
	        
	    end
	
	    function r=results2rise()
	        
	        draws=chain.';
	        
	        r=struct();
	        
	        for id=1:N
	            r.pop(id)=struct('x',draws(:,id),'f',nan);
	        end
	        
	        r.stats.accept_ratio=results.accepted;
	        r.stats.funevals=nan; r.stats.time_elapsed=nan;
	        r.SIG=results.qcov; % covariance matrix now hidden
	        r.c=results.adascale^2;
	        
	    end
	
	end

As the code shows, there are three main articulations of the program :

	#. Obtain from RISE and from the user the inputs that DRAM requires
	#. Run DRAM
	#. Transform the output of DRAM so that it can be further processed by RISE

Obtaining the DRAM inputs
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Some samplers can be run directly using the log-posterior kernel, which can be obtain in RISE through the function
:ref:`pull_objective <rise.pull_objective>`.

The DRAM code on the other hand, requires one function to evaluate log-prior and a separate function to evaluate the
log-likelihood.

In RISE, we can get a handle to the function
that evaluates the prior by using :ref:`log_prior_density <rise.log_prior_density>`. We need to remember though that
RISE stores the parameters vertically while DRAM stores them horizontally. Therefore, before evaluating the prior, we
need to transpose the vector of parameters.

.. code-block:: matlab

    prior = @(theta) log_prior_density(model, theta.');

Where "model" is the dsge model object.	

Likewise we can obtain the function computing the likelihood as follows:

.. code-block:: matlab

    likelihood = @(theta) - model.routines.likelihood(theta', model);

The other elements required by the DRAM code are easier to collect. These include :
the starting vector (x0), the number of draws (n), the covariance matrix (qcov),
additional inputs to the likelihood function if any (data).

Running DRAM
~~~~~~~~~~~~~~~~~~~

With the inputs collected, we can then go ahead and run the DRAM sampler.

.. code-block:: matlab

    [results, chain, s2chain] = DRAM code(model, data, param, options);

Transforming the DRAM output
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

If we want to process the output of DRAM through RISE, we need to transform the DRAM output in a form that RISE
expects.

RISE expects each vector to be stored in a structure alongside the value of the posterior kernel. Unfortunately DRAM
does not return the value of the posterior kernel. This is why we just set to "nan" all elements that cannot be
obtained from the DRAM code.

Processing Posterior draws : The mcmc class
=============================================

.. toctree::

   mcmcHelper

.. # include:: mcmcHelper.rst causes problems of "duplicate label"
   

Visualizing posteriors (and priors)
------------------------------------

See :ref:`Visualizing priors and posteriors<visualPriorsPosteriors>`.

Marginal Data Density computation : the mdd class
====================================================

.. toctree::

   marginalDataDensity

.. # include:: marginalDataDensity.rst causes problems of "duplicate label"